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Abstract 

We compute next-to-next-to-leading order QCD corrections to the gluon-induced 
production cross section of Higgs boson pairs in the large top quark mass limit using 
the soft-virtual approximation. In the limit of infinitely-heavy top quark we confirm 
the results in the literature. We add two more expansion terms in the inverse top 
quark mass to the Mt —)• oo result. Since the l/M^ expansion converges poorly, we 
try to improve on it by factorizing the exact leading order cross section. We discuss 
two ways of doing that and conclude that the finite top quark mass effects shift the 
cross section at most by about 10% at next-to-leading order and by about 5% at 


next-to-next-to-leading order. 

PACS numbers: 14.80.Bn, 12.38.Bx, 14.65.Ha 


1 Introduction 


In the coming years, one of the main tasks in particle physics is the understanding of the 
mechanism of the electroweak symmetry breaking. After the experimental determination 
of the Higgs boson mass, the Higgs potential is fully hxed in the Standard Model. How¬ 
ever, it is very important to independently measure the self-coupling of the Higgs boson, 
which can be obtained from studying the production of Higgs boson pairs. Since the 
corresponding cross section is 0(10^) times smaller than the one for single Higgs boson 
production, Higgs boson pair production poses a challenging problem for the LHC, even 
after the luminosity upgrade around 2020. 

There are a number of phenomenological analyses which investigate the possibility to ex¬ 
tract the self coupling from cross section measurements. First studies have been performed 
more than 15 years ago nm. Since the discovery of the Higgs boson there has been an 
increasing interest in this topic and a number of rehned analyses have been performed, 
see, e.g.. Refs. [IHS]. 

Higgs boson pairs can be produced via the fusion of two partons or vector bosons, via the 
radiation off vector bosons, or in association with heavy quarks. Similar to single Higgs 
boson production, the numerically dominant mechanism is gluon fusion although the 
leading order (LO) contribution is loop-suppressed. Due to the larger Yukawa coupling, 
the dominant contribution comes from top quark loops in the Standard Model. For this 
reason we concentrate in this paper on such contributions. 

For the LO order corrections the exact dependence on the top quark mass and the center- 
of-mass energy is known At next-to-leading (NLO) QCD corrections have been 

computed for the hrst time more than 15 years ago mm in the inhnite top quark mass 
limit using an effective theory. Finite top quark mass effects have been investigated in 
Ref. [12] where a systematic expansion in the inverse top quark mass has been applied 
and a quantitative estimate of the quark mass effects has been provided. It has been 
estimated that they do not exceed (9(10%) of the NLO contribution. Finite top quark mass 
effects have also been considered in Ref. [13] where the exact real radiation contribution is 
combined with the effective-theory virtual corrections. As a result, a reduction of about 
— 10% of the cross section is obtained. We will comment in Section |3] on this issue. 

Within the effective theory also next-to-next-to-leading (NNLO) contributions are avail¬ 
able [T4]IT5] . In this context it is interesting to note that the three-loop matching coefficient 
of the effective operator for two Higgs bosons and two, three or four gluons is different 
from the one for single Higgs boson production ng. The results for the virtual correc¬ 
tions obtained in Ref. [H] have been cross-checked in Ref. [16] where the calculation has 
been performed without reference to the effective theory. The resummation of threshold- 
enhanced logarithms to next-to-next-to-leading logarithmic (NNLL) accuracy has been 
performed in Refs. mm- 

The remainder of this paper is organized as follows: In the next Section we review the 
construction of the soft-virtual approximation for the production cross section and discuss 
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two options to factorize the exact LO result from the higher order contributions. We argue 
that a factorization at the level of the differential cross section w.r.t. the Higgs boson pair 
invariant mass leads to more stable results. Afterwards we reconsider in Section [3] the 
top mass corrections at NLO. Virtual NNLO corrections including finite top quark mass 
effects are computed in Section 01 They are used in Section [5] to present phenomenological 
results for Higgs pair production up to NNLO. Section [6] contains our conclusions. 


2 Factorizing the exact LO expression 


We write the perturbative expansion of the partonic cross section for the production of 
Higgs boson pairs in the form 


aij^HH+x{s,p) = 6ig6jgalg>{s,p) + —(Tlj\s,p) 


etc 



( 1 ) 


where = af\p), ^/s is the partonic center-of-mass energy and ij G {gg-iQg^qg-iqq]- 
Since the quark-induced channels are numerically small m we consider in this paper only 
the gg channel. We use the variable 


P 


m 


H 




( 2 ) 


to parametrize the dependence of the cross section on the Higgs boson and top quark 
mass. We renormalize the top quark mass in the on-shell scheme. Furthermore, we set 
the factorization and renormalization scale equal to each other and write p = pr = Pf- 

For later convenience we introduce for agg^nn+x 

agg-,HH+xis, p) = cr^° + + ... , (3) 


and denote the sum of the first two terms on the right-hand side by 

Finite top quark mass effects io gg ^ HH aX NLO have been considered for the first time 
in Ref. [12] . The applied method is based on “reversed unitarity” [19] where, with the help 
of the optical theorem, the imaginary part of forward scattering amplitudes are computed 
to obtain the total cross section. The virtual corrections have also been computed by 
directly considering the gg ^ HH amplitude. Expansion terms up to order of the 

NLO contribution to a{pp HH) have been computed [T2l[20l[2T]. The factorization of 
the exact LO corrections has then been implemented at the partonic level for the total 
cross section using 
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where cr^g^exact contains the exact dependence on s and p. 

In Ref. m, where the NLO resnlt has been compnted for the hrst time in the heavy top 
qnark limit, a different approach has been applied. Actually, the exact LO result has been 
factorized before the integration over the Higgs pair invariant mass. In this approach the 
(total) NLO cross section reads 


a 


( 1 ) 

99 



dcr, 


( 1 ) 




dcr. 


(0) 


dQ2 


( 5 ) 


where ’exact’ and ’exp’ reminds whether exact or (in p) expanded results are used. is 
the invariant mass squared of the Higgs boson pair. Expressed in terms of da/dQ^ it is 
possible to re-write Eq. (j4]) as 
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( 1 ) 

99 
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(0) 

P'P', exact 
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( 0 ) 

39,exp 



( 6 ) 


From the comparison of Eqs. ([5]) and ([6]) one expects that (jS]) leads to better results since 
the differential factorization (DF) in Eq. (|5]) results in a better-behaved integrand, in 
particular for Q 2 > 4M2. 

For the virtual corrections, which are proportional to 6{s — Q‘^), one has immediate access 
to the Q2 dependence and the DF of Eq. ([S]) can be applied0 The real corrections, 
however, are obtained with the help of the optical theorem which directly leads to the 
total cross section and thus Eq. ([5]) can not be used. For the construction of the soft- 
virtual (SV) approximation, which is discussed below, we need in addition to the virtual 
term also the contributions from soft gluon emission. Since the soft contributions are 
proportional to the LO Born cross section Eq. ((5]) can immediately be applied to the SV 
cross section. 


In the following we discuss the construction of the SV approximation for the cross section 
(see also Ref. [ 22 ] )• For simplicity we consider for the following schematic reasoning the 
total cross section a. Note, however, that the arguments also hold for da/dQ^. in a hrst 
step we split a according to 


a 


^virt+ren 


+ a 


real+split 

5 


( 7 ) 


where the two terms on the right-hand side correspond to the virtual corrections (including 
ultra-violet counterterms) and the real corrections (including the contributions from the 
factorization of initial-state singularities). They are individually divergent and only the 
sum is hnite. In the next step we re-organize the terms on the right-hand side such that 
a can be written as 


— Ssv + Dh , 

^In fact, for the virtual corrections Eqs. o and ([5]) are equivalent. 


( 8 ) 
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where “SV” refers to “soft-virtual” and “H” to “hard”. This is achieved by splitting 
^virt+ren ^ divergent and a hnite term and separating (divergent) soft 

and (finite) hard contribution. The soft contribution is combined with (j™t+ren 
Ssv such that 


I I Y^rea,l-|-split 

SV — ^div + ^fin + ^soft ’ 

Sh = (9) 

Ssv and Eh are separately hnite. Sdiv is constructed following Ref. [23]; explicit expres¬ 
sions can be found in Refs. [T6l|22]. Note that the hnite part is constructed solving the 
equation 


a 


virt+ren 


— Efin -|- E 


div ' 


( 10 ) 


for Efin. In this paper jg computed including top quark mass ehects. Mass ehects 

are automatically taken into account in Ediv and since these contributions are 

proportional to the exact LO cross section [22]. Note that Eqs. (CD, (ED, (ED and oa 
also hold for the diherential cross section da/dQ^ and thus a factorization as suggested 
in Eq. ([5D can be performed for the soft-virtual contribution. 

To be more specihc we present the NLO and NNLO diherential version of Eq. (Iiup which 
reads [T61I22] 


da, 


( 1 ) 


df 

dt 


+ 2Re K‘>] 

df L ^ Wt 

^ + 2Re[/W]i^-- 


+ ||JW| +2Re 


Vt(i)v 


2Re 


dgW 

dt 


( 11 ) 


where t = {qi — with qi (gs) being the incoming (outgoing) momentum of a gluon 
(Higgs boson), a*'°^ = a^° and explicit expressions for the operators can be found 
in Ref. [23] . 

We adopt the notation from Ref. [22] and parametrize radiative corrections to the partonic 
cross section via 


with 



a^^zG{z) , 



( 12 ) 


( 13 ) 


and 

G(z) = i(i-j) + gGW(j)+(^)T'^)(j) + ... 
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= Gsy{z) + G}i{z), 


( 14 ) 


where the renormalization scale dependence in the strong coupling constant ag is sup¬ 
pressed. From Eq. flT^ one obtains for the total cross section 


a = 




dza^^{zs)G{z), 


(15) 


with 


(5 = 



(16) 


In the second line of Eq. ffT4)) we split G{z) into soft-virtual and hard contribution. Note 
that in our approach we do not have access to G]i{z). Actually, at NNLO we only have 
Gsv(z) at hand and at NLO only the heavy top expansion of dz a^^(zs) Gh(z) is 
available to us. 

Explicit results for Gsv(^) can be found in Ref. [22] including higher order terms in 
e specifying, however, the renormalization and factorization scale to ^/s ~ \f^■ For 
completeness we provide the NLO and NNLO results for Gsv('2) for generic p in the limit 
e —0. The results read 
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where Cp = 4/3, Ca = 3 are QCD colour factors, n/ = 5 is the number of massless quarks, 
L = log(/i^/s) and 


D_, = S{l-z) , 


Dn>0 — 


log^(l - 

1-z 


(18) 


The quantities and are evaluated for fi = ^/s. They are obtained from da^^gj^/dt 
and in Eq. flTTl) after integration over t. 

In the next Section we apply the formalism described in this Section at NLO and compare 
to the results obtained in Ref. [12]. The numerical effects at NNLO are presented in 
Section [5] using the mass corrections computed in Section 01 


3 Revisiting NLO 

In this Section we restrict ourselves to NLO and compare the results of Ref. [12] with the 
alternative factorization discussed in the previous section. 

To obtain numerical results we employ MSTW2008 parton distribution functions 
(PDFs) [21] and consistently use N^LO PDFs to compute N^LO {k = 0,1, 2) cross sec¬ 
tions. We assume the energy of the LHC to be 14 TeV and adopt the values of the strong 
coupling constants as{Mz) that we use in our computation from the MSTW PDF £t: 

= 0.13939 , = 0.12018 , = 0.11707. (19) 

We renormalize the top quark in the on-shell scheme and use M* = 173.21 GeV [25] . 
For the Higgs boson we use rriH = 125.09 GeV [26]. For numerical results shown in this 
Section the renormalization and factorization scales have been set to 2mH- 

In Fig. [1] we show the results for the partonic cross section as a function of ^/s from 
Fig. 6 of Ref. [12], see solid lines. The splitting of these results into soft-virtual and 
hard contributions is shown as dotted and dashed curves, respectively. For the complete 
result and soft-virtual approximation the infinite-top quark result corresponds to the (red) 
second line from below. The lowest line includes in addition the 1/M/ corrections. The 
curves including higher order mass corrections are above the infinite-top quark result. 
The topmost curves include 1/M/^ terms. The hard contribution shows a quite flat 
behaviour above ^/s ^ 400 GeV and exhibits smaller shifts when including top mass 
corrections. Furthermore, the infinite top mass result corresponds to the topmost curve 
and the lowest curve includes 1/M/^ terms. From Fig. [1] it is evident that the hard 
contribution is numerically much smaller than the soft-virtual one. 

In Fig. [21 we compare the solid curves from Fig. [1] (which are dotted in this plot) with 
the results obtained with the help of DF applied to the SV approximation and adding 
the hard contribution as given by the dashed lines in Fig. [H The relative position of the 
lines and the colour coding is as in Fig. [H For lower values of y/s the two approaches lead 
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Figure 1: Partonic cross section as a function the center-of-mass energy including various 
orders in the inverse top quark mass. The dotted and dashed curves show the breakup of 
the complete result (solid lines) into soft-virtual and hard contribution. 


to comparable results. However, the DF curves have their maximum at lower values of 
^/s and lead to a smaller cross section for larger values of i/s. Furthermore, one observes 
a drastic improvement in the convergence behaviour when including higher order mass 
corrections. In particular, for y/s = 400 GeV the difference between the inhnite top mass 
result and the one including terms amounts to only 0.05 fb to be compared with 

~ 0.25 fb for the dashed curves. 

At this point we want to stress that the splitting between the soft-virtual and hard 
contribution in Eq. (|8]) is arbitrary. In fact, the soft-virtual contribution of G{z), Gsy{z), 
is constructed for the limit z 1 and thus it is possible to replace Gsv(2:) by f{z)GsY{z) 
with /(I) = 1 (see, e.g., discussions in Refs. [27ll2^ ). The hard contribution is modihed 
accordingly such that the sum of Ssv + Fh does not change. One observes that for 
f{z) = z the soft-virtual contribution approximates the partonic NLO contribution to 
the cross section very welH such that at the hadronic level the deviations are below 2%. 
Based on this observation we use f{z) = z for the NNLO cross section where we only 
have the soft-virtual approximation at hand. Furthermore, better results are obtained by 
replacing L in Eq. (ITT)) by L = which is justihed since ^/s zz in the soft 

limit. 


^Note that the corresponding curves are not shown in Fig. [5] 









SV+H 

SV(DF)+H 


Figure 2: Partonic cross section as a function the center-of-mass energy including various 
orders in the inverse top quark mass. The dashed and solid curves correspond to the 
factorization for the total and differential cross section, respectively. The colour coding 
is taken over from Fig. [T] 


Note that in Ref. [15] it has been observed that the soft-virtual approximation constructed 
in Mellin space approximates the full (effective-theory) result with an accuracy of 2%. 

It is interesting to look at the partonic K factor which is defined via 


^NLO 


^LO _j_ ^(jNLO 


( 20 ) 


Results for the two methods to factorize the exact LO term are plotted in Fig. [3] as a 
function of i/s where the dashed curves are already shown in Ref. [T2|. One observes that 
DF leads to a lower K factor and that the spread among the various p orders is smaller. 
Furthermore, it is interesting to note that for DF the top quark pair threshold behaviour 
of the LO term is not washed out in contrast to the dashed curves. It is common to 
both factorization methods that there is a strong raise when approaching the threshold 
for Higgs boson pair production (see also discussion in Ref. ini)- 

Fig. m shows the hadronic cross section an for Higgs boson pair production including 
NLO corrections as a function of which is a technical upper cut on the partonic 

center-of-mass collision energy. It is introduced via 


O’nisH, Scut) 



('^) 6'(Scut - TSh) , 


( 21 ) 
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SV(DF)+H 

SV+H 


Figure 3: Partonic NLO K factor for the factorization performed at the level of the total 
(dashed) and differential (solid) cross section. 


where the luminosity function is given by 


( 22 ) 


fg{x) are the gluon distribution functions in the MS scheme. Note that in the soft limit 
y^Scut is a good approximation to Q^. The various lines in Fig. 0] correspond to the 
inclusion of different orders in p at NLO. For convenience we show on the right end of 
Fig. mthe total cross section for ^/sh = 14 TeV. Note that the approximation used for 
the computation of the p"’ terms is not valid for large values of \/S(,ut (neither is the 
effective-theory result). However, it can be used as an estimate of the mass correction 
terms. Using the spread as an estimate for the uncertainty we conclude that a hnite top 
mass induces a ±10% uncertainty on top of the inhnite top quark mass result. 

The lower panel of Fig. |4] shows the hadronic K factor which is obtaind from Eq. (I20ll 
by replacing a by an and using NLO PDFs in the numerator and LO PDFs in the 
denominator. j-^igeg close to threshold, however, for ^/Scut ^ 500 GeV one observes 

a flat behaviour of 1.6 (for p = 2m h)- 

Top quark mass effects to double Higgs boson production have also been considered in 
Ref. [13]. In the approximation used in that reference the real corrections are treated 
exactly, however, the inhnite top quark mass approximation is used for the virtual cor¬ 
rections. A decrease of the cross section by about 10% due to hnite top quark mass is 
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Figure 4: NLO hadronic cross section and K factor as a function of 


reported. 

In Fig. [5] we split0 the partonic results for the solid lines of Fig. [1] (which corresponds to 
the factorization according to Eq. (jl])) into virtual corrections (including the ultra-violet 
counterterms; dotted lines) and the real-radiation part which include the contributions 
from the splitting functions (dashed lines). 

We observe that for ^/s ^ 400 GeV, the region where our approximation is valid, the 
top quark mass corrections to the real radiation part (upper dashed curves) reduce the 
inhnite top result by about 10%, in agreement with the observations of Ref. [13]. On the 
other hand, the top mass effects to the virtual contribution leads to a positive shift as 
compared to the effective theory result. Summing real and virtual corrections leads to an 
overall positive effect from top mass corrections as can be seen by the solid curves, see 
also Fig. 121 Note that up to ^/s ~ 400 GeV top mass corrections are dominated by the 
l/Mf terms. 

^Note the the individual terms contain 1/e poles which cancel in the sum. In Fig. [5] the finite contri¬ 
butions are plotted. 
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Figure 5: Splitting of partonic cross section (solid lines) into real (uppser dashed lines) 
and virtual (lower dotted lines) contributions (see text for details). 


4 Top quark mass corrections at NNLO 

In this Section we compute the virtual corrections to the NNLO cross section including 
top quark mass corrections. Afterwards we construct Sgn as described in Section [2] and use 
Eq. (ITTD to evaluate the partonic and hadronic cross sections including 1/Mf correction 
terms. 

4.1 Calculation 

We have applied two methods to compute the virtual corrections. In the hrst one we 
consider the amplitudes for gg —)■ HH up to three-loop order and in the second one 
the forward scattering amplitude is considered, which, after taking the imaginary part, 
directly leads to the total cross section. In both cases we perform an asymptotic expansion 
for large top quark mass. The hrst approach has the advantage that it is straightforward 
to introduce Higgs boson decays whereas the second approach can immediately be applied 
to real corrections. 
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Figure 6: One-, two- and three-loop Feynman diagrams contributing to the process gg —)■ 
HH. Solid lines refer to top quarks, curly lines to gluons and dashed lines to Higgs bosons. 


4.1.1 Amplitude gg —> HH 


NNLO calculations require corrections up to three-loop order to the process gg —)■ HH. 
Typical contributing Feynman diagrams at LO, NLO and NNLO are shown in Fig.|6l They 
are generated with the help of qgraf [22]. Note that in this approach no contributions 
with external ghosts have to be considered since we project on physical states. The 
transformation to FORM [30| code is done with the program q2e [3n[32] and the asymptotic 
expansion for large top quark masses is realized with the help of exp [311 [32]. After 
expanding the identihed hard subgraphs in the small quantities one arrives at one-scale 
vacuum integrals up to three loops and massless one- and two-loop four-point diagrams 
with two massless and two massive external momenta. The vacuum integrals are computed 
using MATAD [33]. In the case of the four-point integrals we apply FIRE [3ll[35] to express 
them as lin ear combinations of master integrals. The latter are shown in Fig. [7] where we 
have qf = q 2 ~ 0^ qi = qt = '^‘h- Analytic results for all master integrals can be 
found in the literature (see, e.g.. Refs. [52H3S])- In this paper we only show results for the 
triangle graph in the second line of Fig. [7] since for our purpose the representations given 
in Refs. [MIEE] are less suited. We use instead 


7i(4) = ) <( Go(-l; a;)Go(0; y) - G'o(0; y)Go{-l/y, x) + 2Go{-l, 0; x) 


H 


- 2Go{-l/y,0;x) -F e 


- 2i7iGo{-l, 0; x) - Go{-l/y; x)Go{0, 0; y) 
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+ Go(-l; 2 :) - i7rGo(0; y) + Go(0, 0; y)^ + 2mGo{-\ly, 0; x) 

+ Go(0; y) {i'nGo{-\ly\ x) + Go(-l, -1; x) + 2Go(-l, 0; x) 

- Go{-l,-l/y,x) + Goi-l/y, -l]x) - 2Go{-l/y,^]x) - Go{-l/y, -l/y]x)^ 

+ 2Go(-l, -1, 0; x) + 4Go(-l, 0, 0; x) - 2Go(-l, -l/y, 0; x) + 2Go{-l/y, -1, 0; x) 


- 4Go(-l/|/,0,0;a;) 


2Go(-l/|/, -l/2/,0;x) 



(23) 


with 



_ 1 + \/5 

“ 1 - 
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(24) 


Go{.{wi}; z) are Goncharov Polylogarithms [31] with weight {tCj} and argument 2 ; defined 
through 


Go(«;i,m;2, ...,Wn]x) 
with Wi,x G C and 



dt- - Go{w 2, ...,Wn,t), 

t — Wi 


(25) 


Go(0^, x') 


n\ 


log” X . 


(26) 


The functions of the e° term in Eq. (125]) can be expressed in terms of logarithms and 
dilogarithms via 


Go{0-,y) = \og{y) , 

Go{-l;x) = log(l + a:) , 

Go{-l/y, x) = log(l + xy) , 

Go(-l, 0; x) = log(x) log(l + x) + Li 2 (-a:) , 

Go{-l/y, 0; x) = log(a;) log(l + xy) + U 2 {-xy). (27) 

We have cross checked the numerical result for the e° and terms of /i(4) in Eq. fl23|) 
against FIESTA |1D]. 

Using this method we have computed terms up to order 1/M/^ at NLO [T2ll20| [2T] and 
terms up to order 1/Mf at NNLO. As an important check we have computed the 1/Mf 
corrections for general QCD gauge parameter which drops out in the final expression. 

From the calculation of gg —)■ HH one obtains in a first step results for da/dt. Integration 
over phase space then leads to da/dQ"^. For the results in Section [5] this integration is 
performed numerically. 
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Figure 7: One- and two-loop master integrals needed after applying asymptotic expansion 
to the amplitude gg —?■ HH. All internal lines are massless, qi = 02 = 0; = qt = 


4.1.2 Amplitude gg — )■ gg 

The second method is based on the use of the optical theorem in analogy to the NLO 
calculation performed in Ref. [12]. This method serves as an important cross check. In 
the following we provide some of the technical details 

• The amplitudes for gg gg are generated with the help of qgraf [29] . 

• In a hrst step about 17 million diagrams are generated. However, most of them do 
not contain a cut through exactly two Higgs bosons. For this reason we post-process 
the qgraf output and hlter [2T|ITT] the amplitudes describing the virtual corrections 
to gg ^ HH. Typical Feynman diagrams are shown in Fig. [HI 

• FORM [2011121 code is then generated by passing the output via q2e [3I1132], which 
transforms Feynman diagrams into Feynman amplitudes, to exp [211122] . 

• Our in-house FORM code applies projectors {—g^Mu) for each pair of external gluons 
which includes also non-physical degrees of freedom in the sum. Thus also contri¬ 
butions with ghosts in the initial state have to be considered. Note that this is in 
contrast to single Higgs boson production which has no virtual contributions with 
ghosts in the initial state. 
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Figure 8: LO, NLO and NNLO Feynman diagrams needed for the forward scattering 
amplitude gg ^ gg. Solid lines refer to top quarks, curly lines to gluons and dashed lines 
to Higgs bosons. At NNLO only virtual contributions are shown. Wavy lines denote the 
cuts. 



Figure 9: Resulting Feynman diagrams after shrinking the top quark loops to a point 
according to the rule of asymptotic expansion. The blob represents one-loop vacuum 
integrals. 


The application of the asymptotic expansion for large top quark mass leads to a factor¬ 
ization of the hve-loop integrals into the following contributions: 
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Figure 10: Phase space master integrals occurring in the amplitude for gg —)■ gg. The 
hrst line contains LO and NLO integrals. The integrals in line two and three are needed 
for the NNLO virtual corrections. Single and double lines represent massless and Higgs 
propagators, respectively. Double lines with gray-shaded interspace (last two diagrams in 
hrst row) correspond to Higgs boson propagators which shall not be cut. A cross marks 
an inverse propagator. 


1. four one-loop vacuum integrals times one-loop integrals which are already present 
at LO, see Fig. |9](a); 

2. three one-loop vacuum integrals times two-loop integrals, see Fig. |9](b); 

3. two one-loop vacuum integrals times three-loop integrals, see Fig. |9](c); 

4. two-loop times one-loop vacuum integrals times two-loop integrals; 

5. three-loop times one-loop vacuum integrals times one-loop integrals. 

The mass scale in the vacuum integrals is given by the top quark. They are again computed 
with the help of MATAD [33] . For the remaining integrals, which depend on <5, we use the in- 
house programs rows [H] and TopoID [^ITT] to perform the reduction to master integrals. 
The latter are depicted in Fig. HU] 

All three-loop master integrals factorize into a two-loop form factor contribution and the 
one-loop master integral of the LO calculation. From the latter only the imaginary part 
is needed which is well known. The results for the two-loop form factor integrals can, e.g., 
be found in Refs. [43] . 
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The two-loop master integrals are more involved. A numerical calculation would probably 
be possible, however, we follow the approach outlined in Ref. [12] for the NLO master 
integrals and perform an expansion in 5 = 1 — Amjj/s [cf. Eq. flTHl) ]. All integrals contain 
a massless sub-loop for which analytic expressions are known. The massless two-point 
function can be expressed in terms of T functions and results for the triangle with two 
massless external legs and the (crossed) box can be found in Ref. |37| . Analytic expressions 
for the triangle diagram with squared external momenta s, rn\j, are given in Eq. fl23|) . 
After expanding in 5 the remaining phase-space integration can be performed analytically. 
We have computed expansion terms up to order 5^° and found agreement with the results 
obtained in the previous subsection which for this purpose have also been integrated 
analytically after performing the expansion in 5. 

The approach based on the optical theorem requires special care in the treatment of the 
imaginary parts originating from the two-loop form factor diagrams. Such contributions 
either correspond to or In the former case the 

two-loop integrals originate from the product of two one-loop contributions containing 
factors (—l-(-fO)'^ and ( —1—iO)^, respectively, which hnally leads to (— 1 -|- 20 )'^(—1—iO)'^ = 1. 
In the other case one has (—1(—1 — = 1 — -|-(!2(e^). The corresponding 

discussion for single Higgs production can be found in Ref. [Hj. 

Note that the approach based on the gg ^ HH amplitudes leads to simpler intermediate 
expressions. Thus, it is possible to allow for a general QCD gauge parameter ^ when 
computing the 1/M| terms. Furthermore, also the l/M^ corrections could be evaluated 
(for = 0) whereas in the optical theorem approach only l/M^^ terms could be computed 
in Feynman gauge. However, let us mention that this approach can be used in a straight¬ 
forward way to compute NNLO top quark mass effects to the hard contribution of the 
total cross section whereas in the approach of the previous subsection this is less obvious. 

To conclude this Section let us summarize our procedure to obtain the SV corrections at 
NNLO. We compute virtual corrections to gg ^ HH including 1/Mf corrections. Note 
that we have dcr/dQ^|virt ~ ~ s)- Using Eq. ffTUD we construct crgn which enters Gsv 

in Eq. (1T7I) . The differential and total cross section is then obtained with the help of 
Eqs. (fT2|) and ([T5|). 


5 Improving NNLO 

In this Section we discuss the effect of the l/M^ and l/Mf terms on the NNLO cross 
section for the production of Higgs boson pairs. The cross section in the infinite top 
quark mass limit has been computed in Ref. [Tl|. In Ref. (TB] the three-loop matching 
coefficient has been added, completing the NNLO prediction, and the virtual corrections 
from Ref. [T3| have been cross checked. 

The results for the virtual corrections computed in Section 0] are inserted in the formalism 
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Figure 11: Comparison of the LO, NLO and NNLO contributions to the partonic cross 
section. At LO the exact result is shown and at NLO and NNLO the hrst three terms in 
the large-expansion are shown. For all curves the NNLO-value for is used. For the 
renormalization and factorization scale we use p = 2m h- 


described in Section [2] to construct the quantity Ufin which enters Eq. flTTI) . The result 
for the partonic cross section is shown in Fig. [TT]as a function of the partonic center-of- 
mass energy where the exact LO result is compared with NLO and NNLO. At NLO and 
NNLO three terms in the mass expansion are shown0 Furthermore, at NNLO the SV 
approximation is shown for f{z) = z (cf. discussion in Section [3]). Note that the NNLO 
curves peak for smaller values of y/s than at NLO and LO. As far as the top quark mass 
corrections are concerned the same pattern is observed as at NLO: the correction term of 
order p decreases the inhnite top quark mass result which is overcompensated by the 
term resulting in a small positive correction. 

In Fig. [12] we compare the NNLO-SV contribution for two different scales, p = 2mH and 
h = \/^- This plot furthermore shows the effect of f{z) = 1 and f{z) = z. Note that 
the choice f{z) = z, which we expect to better approximate the complete result, leads to 
an increase of the cross section. 

The hadronic cross section as a function of is shown in Fig. [13] for p = 2m h- At LO 
the exact result is used and the NLO curve is based on the results. Using instead the 
p® terms leads to an upwards shift of about 5% as can be seen in Fig. ]1] At NNLO three 
curves are shown which include terms up to p°, p^ and p^. As at NLO one observes good 
convergence up to ~ 400 GeV. For higher values of the p^ curve is below and 

^In this plot we only include mass corrections up to order at NLO to have a direct comparison with 
NNLO. 
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Figure 12: NNLO partonic cross section for different scales and for f{z) = 1 (dotted) and 
f{z) = z (solid). Only the result is shown. 



250 350 450 550 650 750 oo 

[GeV] 

Figure 13: Hadronic LO, NLO and NNLO-SV (with f{z) = z) cross sections as a function 
of y'Scut- For their evaluation the respective value of as is used. At LO the exact result 
and at NLO only the term is shown. At NNLO the p^, and p^ results are plotted. 
The results in the right panel with “cxo” at the bottom correspond to the prediction of 
the total cross section. For this plot p = 2mH has been used. 
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Figure 14: Hadronic NLO (dotted) and NNLO (solid) K factor as a function of y'Scut- 


the curve above the inhnite top quark mass result leading to a ±2.5% effect for the 
total cross section on the rightmost part of the plot. To be conservative, we thus estimate 
that the NNLO top quark mass effects lead to an uncertainty of ±5%. Note that the 
NNLO corrections amount to about 20% of the LO result. 

Fig. im shows the hadronic K factor at NNLO which is dehned by 

I A^NLO I r^NNLO\ | 

I^NNLO _ jiNNLOpdfs 

- “XOI • 

H IlO pdfs 


as a function of \/Scuf comparison also the NLO result from Fig. |4]is shown as dotted 
curve using the —)■ oo result. One observes that the various p orders lead to similar 

results for Furthermore, there is a strong raise close to threshold which is due to 

the steeper behaviour of the NNLO cross section as can be seen in the inlay of Fig. [HI 
For higher values of ^Scut, in particular for the total cross section, approaches 

1.7- 1.8. 

Results for the total cross section at LO, NLO and NNLO are summarized in Table [T] for 
p = 2mH- At NLO only terms are included in the analysis whereas at NNLO p^, p^ and 
p^ terms are considered. In this way we can estimate the top mass effects of the NNLO 
term. Besides the cross section also the K factor is shown. At NNLO we use f{z) = z 
and we furthermore show the relative deviation due to 1/M|"' terms dehned through 


^NNLO 




a 


NNLO I 
H I 


(29) 
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an [fb] 

^(N)NLO 

^NNLO 

LO 

22.7 

— 

— 

LO±NLO|p0 

36.4 

1.60 

— 

LO±NLO p0±NNLO|p0 

39.7 

1.75 

0 

LO±NLO poTNNLO pi 

38.7 

1.70 

-2.5 

LO±NLO p0±NNLO|p2 

40.5 

1.78 

±2.0 


Table 1: Total hadronic cross section at LO, NLO and NNLO-SV including top quark 
mass effects using /i = 2m h and f{z) = z. 

The two known mass correction terms lead to a change of the cross section by about 
±2%. Assuming a similar pattern as at NLO we thus estimate that NNLO top quark 
mass corrections change the effective-theory result by at most ±5%. 

6 Conclusions 

We compute NLO and NNLO corrections to double Higgs boson production in gluon 
fusion beyond the effective-theory approach. The starting point of the calculation are 
full-theory Feynman diagrams. We perform an asymptotic expansion in the limit where 
the top quark mass is large and compute at NNLO three terms in the l/M^ expansion for 
the virtual corrections. They are used to construct a soft-virtual approximation for the 
production cross section. In the limit —)■ cx) the effective-theory result of Ref. [H] is 

conhrmed and and 1/Mf terms are added. 

The main result of this paper is illustrated in Fig. [13] where the hadronic cross section is 
shown as a function of y'Scut (a technical cut on the partonic center-of-mass energy). The 
curves including 1/Mt corrections deviate from the infinite mass result only by a few per 
cent which leads us to the estimate that the effective-theory result is accurate to ±5%. 
Analog results for the mass corrections at NLO are shown in Fig. |4| which constitutes an 
update of Ref. [12]. Here we estimate the uncertainty to ±10%. 

We want to stress that the results obtained in this paper provide excellent approximations 
for small values of the partonic center-of-mass energy, say below y/s ~ 400 GeV. Although 
in this region the cross section is small it is of interest since there the cross section has a 
characteristic behaviour. Furthermore, it is possible to use our result in this region as a 
benchmark for future calculations taking into account the exact dependence on Mf. 

The methods described in Section ID can also be used to compute top mass corrections 
to the real radiation part. However, the simplifications used in Ref. [15] where results 
have been obtained for Mt —)■ cx do not apply once finite mass effects are considered. 
The calculation is much more challenging since significantly more Feynman diagrams 
contribute and more complicated master integrals have to be computed. 

In this paper for the first time the effect of a finite top quark mass has been examined 
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for the NNLO cross section for double Higgs boson production. Whereas at NLO an 
exact calculation is within reach this is certainly not the case at NNLO. Thus our results 
become particular important once our NLO approximations are compared to an exact 
calculation which increases the confidence in the uncertainty estimate. Furthermore, one 
probably can obtain a prescription to tune the approximation procedure and hence reduce 
the uncertainty at NNLO. 
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